Correlation density matrix: an unbiased analysis of exact diagonalizations 
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Given the ground state wavefunction for an interacting lattice model, we define a "correlation density matrix" 
(CDM) for two disjoint, separated clusters A and B, to be the density matrix of their union, minus the direct 
product of their respective density matrices. The CDM can be decomposed systematically by a numerical 
singular value decomposition, to provide a systematic and unbiased way to identify the operator(s) dominating 
the correlations, even unexpected ones. 
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The ground state of a strongly-interacting, quantum- 
mechanical lattice model (with spin, boson, or fermion de- 
grees of freedom) is characterized by long-range order, power- 
law correlations, or the lack of these. When such a system is 
studied numerically, it may be unclear a priori what kind of 
correlation will be dominant - especially in cases where exotic 
order or disorder are possible, such as the doped square-lattice 
Hubbard model, or (better) the highly frustrated s = 1/2 
Kagome antiferromagnet; in the latter system spin-spin, spin- 
Peierls, spin-nematic, or chiral order parameters were all seri- 
ous candidates OJ]. Before computing the ground state corre- 
lations, one must first guess which operators are important - a 
choice which is necessarily biased by one's prior knowledge 
or preconceptions, and is problematic for hidden or exotic or- 
ders. 

In contrast, approaches based on the density matrix (DM) 
of a cluster of several sites are unbiased - apart from specifi- 
cation of that cluster - since the DM specifies the expectation 
of every operator local to the cluster - including the "key op- 
erator^)" meaning those having long range order (i.e. order 
parameter) or having strong correlations. For exact diagonal- 
izations (ED) of interacting systems, the DM was used as a 
diagnostic to compare different system sizes |2[] or truncations 
of the Hilbert space (U . 

Here we propose a new application of the density matrix as 
a way to uncover correlations/orders from numerics without 
requiring any foreknowledge of what kinds to expect. Con- 
sider two small disjoint clusters A and B (identical apart from 
a translation), either cluster having a Fock-Hilbert space of 
dimension D. Let p AB be the many-body density matrix for 
the disconnected "supercluster" A U B, constructed from the 
whole system's ground state wavefunction by tracing out all 
other sites, with p A and p B similarly defined. Then we define 
the correlation density matrix (CDM) to be 



f EE p AB 



P A 



(1) 



If there were no correlations between clusters A and B, then 

pAB _ pA ^pB an( j pC _ q 

The CDM defined in (Q3 contains all possible inter-cluster 
correlations J3]. Write the ("connected") correlation of the 



fluctuations of any two operators as (PQ) C = (PQ) — 
(P)(Q); then if P(A) and Q(B) act on clusters A and B, 



(P(A)Q(B)) c = Tr[p c P(A)Q(B)}. 



(2) 



Index relabeling and the operator singular-value decom- 
position — The key notion underlying our processing of the 
CDM is, given the D x D matrix representing an operator 
on a cluster's D dimensional Hilbert space, to rewrite it as 
an D 2 -component vector of complex numbers using fused in- 
dices H] (a', a) <-> a(a\ a), (&', b) <-> f3(b', b). Say that p c is 
known in terms of the product states \a')\b') and \a) \b) of the 
occupation-number basis on the clusters [6]. Then 

P C = E P%, ab W)W){a\{b\^Y, C ^9 a h (3) 



where p% . = C c 



(a' ,a),P(be,b)- 



a/3 



Here g 



\a')(a\ and 



hp ee \b')(b\ are bases for the respective clusters A and B, 
manifestly orthonormal in terms of the Frobenius norm 



\Hf 



J2\Pa',a\ 2 =Tr(P^P 



(4) 



for any operator P, and the Frobenius inner product 

(P, Q) F ee £ P:^ a Q al . a = Tr (ptQ) . (5) 



(In the fused-index notation, Eqs. (0]i) and (O take on the usual 
form of a vector norm and vector inner product.) 

Next a numerical singular value decomposition can be 
made of C a p as a matrix of complex numbers: 



C a /3 = E GvUvaVv 



■P 



(6) 



where U and V are unitary matrices, and {er„ : v = 
1, . . . , D 2 } are the singular values. [Eq. © can also be writ- 
ten in the matrix form C = U T YA/, where S = diag({<7„}).] 
Substituting (|6]l into (O, we obtain the operator singular- 
value decomposition, 
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D 



P C = ^2^X V (A)%(B) 



(7) 



FIG. 1 : Model: spinless fermions, with hardcore excluding nearest- 
neighbors, on a ladder, with longitudinal hopping tn = 1, trans- 
verse hopping t±, and correlated hopping t' . The correlation density 
matrix involves two clusters, each of 2 x 2 sites, with their centers 
(marked +) separated by r. This ladder has length L = 8, with peri- 
odic boundary conditions as indicated by the + at right edge. 



This (simple but powerful) expression is the key formula 
of our paper. Each term represents the correlated quantum 
fluctuations of Frobenius-orthonormalized basis operators 0], 
X v = J2 a U va g a on cluster A and Y u = J2p V v php, on clus- 
ter B. 

Recalling @, we can rewrite any correlation 

{P(A)Q(B)) c = J2a v (Xt,P(A)) F (Yj,Q(B)) F (8) 

in terms of Frobenius inner products (3). In particular, 
(l„(A)ty T (B)t) c = a v 6 UT , Thus {X V (A)*} and {Y v (B)i} 
are the natural bases into which operators P(A) and Q(B) 
should be decomposed. Each \<j v \ is a normalized measure 
of the strength of the corresponding inter-cluster ground state 
correlation. By convention, we order the singular values 
0"i > 02 — ' ' ' — a D 2 > 0. This ordering gives a means 
of approximating p c by retaining just the first few terms in 
the expansion ((T). 

Observe that Hp* 7 )! 2 = °v| 2 ls a basis-invariant mea- 
sure of the total correlations between A and B. Since [8] 



\P C \ 



= \\P AB \\1 



\P A \\ 2 F 



(9) 



it follows that 1 1 p c \\ 2 F < 1 — 1/D 2 » 1, which gives a standard 
of comparison for numerically obtained a v 's. 

The CDM typically inherits various symmetries from the 
input wavefunction (ultimately from the Hamiltonian), such 
as spin-rotations, lattice rotations/reflections, or fermion num- 
ber conservation [|9|]. The matrix C a p breaks up into 
symmetry-labeled blocks, which (as with diagonalization) can 
be singular-value-decomposed independently. Each term in 
the expansion (Q is thus assigned to a sector according to the 
quantum numbers carried by X v and Y v , and each sector is 
interpreted as representing a different kind of orrelation. 

A convenient test bed to study CDM properties is a non- 
interacting system (including BCS states) for which density 
matrices can be calculated exactly, iflOn . We analytically 
checked the CDM and its operator SVD for a free Fermi sea 
in one dimension (Ref. [Ill chapters 5 and 6), finding the ex- 
pected FL correlations with an r -1 / 2 envelope and CDW cor- 
relations with an r -2 envelope. 

Ladder model: limiting regimes and operator classes — 
We now test the CDM method on a toy system (Fig. [TJ in 
which spinless fermions hop on a two-leg ladder of length L; 




FIG. 2: Each plot shows (on a log scale) the magnitude of the largest 
singular value for each symmetry sector of operators. The symme- 
tries are labeled "CDW" for number operator (or any combination 
c\cj in the same cluster); "FL" for single creation/annihilation (i.e. 
the correlation function is a 2-point Green's function); "SC" for su- 
perconducting (combination cjct in same cluster). The symmetry 
label ± denotes even/oddness under exchanging the legs of the lad- 
der. In every case, there are 4 particles on a ladder of length L = 8, 
and twist boundary condition averaging was used. (a). No-passing 
ladder with t± = 0.1, if = 0; (b). Rung-fermion case (each fermion 
delocalized on a rung) with t± = 100, t' — 0; SC singular val- 
ues do not appear since they are ~ 10 -15 . (c). Boson pair state: 
t± =0,t' = 100. 



they are forbidden to occupy adjacent sites (i.e., the nearest- 
neighbor repulsion is V — oo). Three kinds of hopping 
amplitudes appear: tn = 1 along legs, t± along rungs, 
and if a "correlated hop" conditioned on a second fermion, 



-t'(c]a 



z^Ci + qcj)n/s; here i, j are two steps apart on the same 
leg, and fik is the number operator for the site between i and 
j on the opposite leg (which would block the tn hops). 
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The phase diagram (see Ref. Fig. 8.1) may be un- 
derstood through the three limiting cases in which one hop- 
ping dominates, (a) tit dominant ("no-passing" limit): the 
leg index is a conserved flavor; the model reduces to a free 
fermion chain (with fermions on alternate legs) (b) t± dom- 
inant ("rung-fermion" limit): each fermion delocalizes on 
a rung, so at low energy the model maps to reduces to a 
fermion chain with nearest neighbors excluded; (c) t' dom- 
inant ("paired" limit): fermions bind into effective (p-wave) 
boson pairs (in one dimension, with nearest neighbors ex- 
cluded). Regime (c) must be dominated by superconductivity 
at large length scales. 

Each of the three limiting cases maps nontrivially to free 
fermions. Elsewhere [ 12] we derived from these maps a semi- 
analytic method ("intervening particle expansion") to calcu- 
late various correlation functions; the results of Ref. [13 have 
illuminated the present calculation. The asymptotic behaviors 
(as expected) are that of a Luttinger liquid: power-law decays, 
with possibilities of commensurate locking when the filling is 
a rational fraction. 

We performed exploratory exact diagonalizations using pe- 
riodic boundary conditions, with four fermions on a ladder 
of length L = 8, the smallest (nontrivial) case at 1/4 filling. 
(This is the most interesting filling - and the hardest, since 
the Hilbert space is largest at filling 1/4: see Ref. 13(b), ap- 
pendix.) The largest block matrix for a sector is 27 x 27. (As 
in our earlier ED studies on the square lattice QdH], the spin- 
lessness and the neighbor exclusion greatly limit the Hilbert 
space compared to e.g. a Hubbard system of the same dimen- 
sions.) To minimize finite-size effects on the density matrices, 
it was necessary to use phase-twist boundary conditions lfl4ll 
(i.e. to thread flux through the "ring" of sites) and average 
over 21 distinct phase angles. (See Ref. 2 and Sec. 8.2.4 of 
Ref.HD). 

Each of our two clusters is 2 x 2 (two adjacent rungs) as 
shown in Fig.[T| the smallest cluster that can capture supercon- 
ducting correlations; each cluster's Hilbert space has dimen- 
sion D = 7. The operators {X Ul Y„}, emerging from the op- 
erator singular-value decomposition, are classified into three 
main categories, according to the fermion number change AF 
they carry: (i) CDW (charge-density-wave-like), AF = 0, 
e.g. the number operator hi on site i H 1 511 z (ii) FL (Fermi- 
like), AF — ±1, e.g. the operator c\ on a site. The two- 
point Greens function, the dominant long-range correlation in 
a Fermi sea, belongs with this operator sector, (iii) SC (su- 
perconducting), AF — ±2; such operators are the order pa- 
rameters for superconductivity. In addition, each operator can 
be even or odd under exchange of the ladder's legs, which we 
denote by appending "+" or "— ". 

Generically, the basis operators {X V ,Y V } do not take the 
minimal form one would adopt in defining a correlation func- 
tion (even in the free fermion case). Instead, complicated 
terms are admixed iflrSll . For example, the dominant opera- 
tor in the FL sector not only has single creation operators cj, 
but terms c\hj. 

Numerical results and conclusions — Fig. |2] presents the 
numerical singular values for the CDM in the three limits; the 



TABLE I: Correlation behaviors in limiting-case models Row la- 
bels (a, b, c) correspond to the panels in Fig. [2] Columns "Sim" 
summarize behaviors inferrable from Fig. [2] "large", "medium", or 
"small" indicate singular values roughly constant with r, i.e. possible 
long-range order (values over 10~\ 1CT 2 , or 10~ 3 , respectively). 
Singular values decaying with r are labeled "d(fast)" or "d(slow)". 
Columns are labeled by the symmetry sectors as in Fig. [2] For com- 
parison, the columns "Th" are from semi-analytic computations of 
Ref. u% exp = exponential decay, LRO = long range order. For the 
pairing limit (c), the FL correlation exponent varies with filling n, 
with a(n = 1/4) « 1.1. 
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d(slow) exp 
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decay behaviors of the different correlations are summarized 
in Table|I] where they are compared with our knowledge from 
the intervening particle expansion fl^tl . Due to the limited 
system sizes for ED, the CDM analysis cannot determine the 
dominant kind of correlation at large distances. That is prac- 
tically impossible for Luttinger liquids in any case: for the 
hardcore boson chain (related to our models) the asymptotic 
(superfluid) correlations may dominate only after 50-100 sites 
B17I1 . TableUshows there is a general correspondence between 
the decay rate of known correlations and that of the singular 
values; the degree of correlation in Fig.|2]tends to be overesti- 
mated due to the very small range of r. 

The rung-fermion case (b) at filling 1 /4 breaks translational 
symmetry, with period-2 long-range order. Examination of 
Fig. |2] (b) indeed shows the corresponding contrast with the 
other two cases: the singular value for the order-parameter 
operator (CDW+) is non-decaying and saturates the bound 
a = 1/2, whereas other kinds of singular values are orders of 
magnitude smaller. 

In the boson-pair case (c), as t' grows large (the boson-pair 
limit), a crossover is expected to asymptotic superconducting 
(SC) correlations; but Fig. [2jc) shows that CDW correlations 
still dominate at all accessible distances, similar to hardcore 
bosons 1 17]. A partial success the CDM analysis is that the 
SC singular values are visibly larger than in the other cases, 
competitive with FL correlations; absent any other knowledge 
of this system, the SC order parameter would be flagged for 
further study (e.g. analytic, or by quantum Monte Carlo). 

In all three cases, most correlations decay generically [12] 
as C(r) ~ cos(2m/ci?r + (S)/^! 2 , where 2mfcj? is an even 
multiple of the Fermi wavevector and x is some correlation 
exponent. Over a small range of r, the with oscillations with 
r obscure the asymptotic r dependence of the singular val- 
ues. We conjecture each such correlation is associated with 
a pair of singular values, oscillating 90° out of phase inside 

/ 1 /2 

the same envelope. Ideally, then, one should plot \J2' lJ <r^] , 
where " runs over just one symmetry sector, to obtain a 
mono tonic decay as l/lrj^. In practice, for reasons we do not 
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understand, this gave little or no improvement. 

To conclude, we have introduced a new tool for analyzing 
exact-diagonalization ground states, using the density matrix 
of a pair of clusters to extract all their correlations in an un- 
biased fashion. Furthermore, via a singular-value decompo- 
sition, the kind of operator dominating the correlations could 
be identified, using (0. There are two regimes where asymp- 
totic decays are not at issue and the correlation density matrix 
based on exact diagonalization should be effective. First, for 
systems believed to have negligible correlations beyond the 
nearest neighbor - e.g. quantum spin liquids in highly frus- 
trated antiferromagnets fljj] - the CDM is the foolproof way to 
confirm the absence of any correlations. Secondly, in systems 
having long-range order [such as our case (b)], the CDM de- 
tects the symmetry breaking. On the other hand, critical states 
[such as the Luttinger liquids of our cases (a) and (c), above] 
are the least promising systems for study by CDM, so long as 
the system sizes are limited by dependence on ED. But if the 



CDM and density-matrix renormalization group methods are 
married 11811 . the asymptotic scaling may become accessible 
for one-dimensional systems. 

Another unbiased method has been proposed to discover 
the symmetry breaking operator from ED using the density 
matrix [19]. It differs from the CDM in two ways: (i) it 
is based on the DM of just one cluster; (ii) it requires not 
only the ground state's wavefunction, but that of several low- 
lying eigenstates which are conjectured to be linear combina- 
tions of symmetry broken states (and degenerate in the ther- 
modynamic limit). That method is meant only for cases of 
long-range order, whereas in principle the CDM identifies the 
strongest correlations even in disordered phases. 
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